We will fit a multi-level non-linear (on the parameters) regression model to the data from the sample. The regression takes the form:

\[ \text{sooner_choice}=\frac{\text{endowment_later}\cdot(b^{t_0}d^k\text{pratio})^{\frac{1}{a-1}}}{1+\text{pratio}\cdot(b^{t_0}d^k\text{pratio})^{\frac{1}{a-1}}} \] where \(\text{endowment_later}=20\).

In R, this looks like the following formula for the population-level parameters (i.e., the fixed effects):

sooner_choice ~ ( endowment_later * (b^t0 * d^k * pratio)^(1/(a-1)) ) / ( 1+pratio * (b^t0 * d^k * pratio)^(1/(a-1)) )

We start by importing the data and ensuring that it contains all of the information for each trial we need. For example, we need to add columns that specify the delay period for each trial.

mctb_data_ <- utils::read.table(file.path(data_dir_base, 'mCTB_kit/data.txt'), 
                         header = TRUE, 
                         na.strings = '.')

mctb_data_l <- stats::reshape(mctb_data_, 
                              varying = paste0('c', 1:24),
                              v.names = 'c',
                              timevar = 'budget_number',
                              times = 1:24, 
                              direction = 'long')

mctb_design <- read.table(file.path(data_dir_base, 'mCTB_kit/instrument_details.txt'),
                     header = TRUE)

mctb_data <- merge(mctb_data_l, mctb_design, 
                   by = 'budget_number', 
                   all.x = TRUE, all.y = TRUE)

head(mctb_data)
##   budget_number subject_id c id endowment_soon endowment_later
## 1             1          1 6  1             19              20
## 2             1          2 6  2             19              20
## 3             1          3 1  3             19              20
## 4             1          4 1  4             19              20
## 5             1          5 1  5             19              20
## 6             1          6 3  6             19              20
##   sooner_date_weeks delay_weeks
## 1                 0           5
## 2                 0           5
## 3                 0           5
## 4                 0           5
## 5                 0           5
## 6                 0           5

We will also take this opportunity to examine missing reponses.

n_missing_responses <- dplyr::summarize(group_by(mctb_data, subject_id),
                                        n_missing = sum(is.na(c)),
                                        p_missing = n_missing/n())
plot(n_missing_responses$subject_id, n_missing_responses$p_missing, 
     type = 'h', xlab = 'Participant ID', ylab = 'Proportion missing responses',
     ylim = c(0,1))

knitr::kable(filter(n_missing_responses, p_missing > 0),
             digits = 2,
             caption = paste0('Proportion missing by participant. Total N with missing data = ',
                              sum(n_missing_responses$p_missing > 0)))
Proportion missing by participant. Total N with missing data = 9
subject_id n_missing p_missing
110 3 0.12
217 1 0.04
270 20 0.83
395 24 1.00
613 20 0.83
624 20 0.83
655 6 0.25
1007 2 0.08
1068 2 0.08

Here, we define some trial-level predictors. t0 defines whether the soonest choice is today or also in the future. We modify k so that it is in terms of days, rather than weeks. We also calculate the price ratio, which is just the ratio of the later reward and sooner reward.

mctb_data$t0 <- as.numeric(mctb_data$sooner_date_weeks == 0)
mctb_data$k <- 7 * mctb_data$delay_weeks
mctb_data$k_w <- mctb_data$delay_weeks
mctb_data$pratio <- mctb_data$endowment_later / mctb_data$endowment_soon

Here we calculate the possible value choices the participant saw for each trial.

for(j in 1:6){
  mctb_data[paste0('soon_', j)] <- (20/as.vector(mctb_data$pratio)) - (j-1) * (20/as.vector(mctb_data$pratio)) / 5
  mctb_data[paste0('late_', j)] <- (j-1) * 4
}

mctb_data$soon_6 <- 0

Finally, based on the particpant’s response, 1-6, we create a column that contains the amount of money they choose to receive at the soonest possible time. If they choose option 6, to receive the full $20 at the later time, they always receive $0 at the sooner time.

#Create a new variable with a value that corresponds to the response option
#chose in `mctb_data$c`
mctb_data$sooner_choice <- apply(
  mctb_data, 1, #for every row in mctb_data
  function(a_row){ 
    #write the name of the column to get the choice value from
    choice_col <- paste0('soon_', a_row['c']) 
    #get the choice value
    sc <- a_row[choice_col] 
    #if it's null, return NA, otherwise return the choice value
    return( ifelse(is.null(sc), NA, sc)) 
  })

Exclusion criteria are specified:

If participants miss greater than or equal to 50% of items on a personality scale, they will be coded as missing a score for that scale and not be included in analyses involving that scale. Due to the nature of the CTB task, participants who miss more than one item per timeframe pair will be coded as missing and excluded from analysis. Finally, participants who did not self-report household income will be excluded from analysis.

As we saw above, some participants are missing more than one response on the CTB task. No one did not report income. Some participants failed to answer >50% of items on a personality scale.

scale_data <- read_csv(file.path(data_dir_base, 'JRPdataOSF.csv'))
## Parsed with column specification:
## cols(
##   .default = col_double(),
##   annualincome = col_character()
## )
## See spec(...) for full column specifications.
scale_data_missing <- select(scale_data, contains('score'), subID) %>%
  filter_at(.vars = vars(-subID), any_vars(is.na(.)))

missing_ids <- full_join(filter(n_missing_responses, n_missing > 0),
                         scale_data_missing,
                         by = c('subject_id' = 'subID'))

mctb_data <- anti_join(mctb_data, missing_ids,
                       by = 'subject_id')

Total number of excluded IDs is 21, or 1.9%.

We use fixed-effects non-linear least-squares as a first pass on the data and to set starting values for the individual-level analyses that will assist the multi-level estimation.

sd_of_soonerchoice <- sd(mctb_data$sooner_choice, na.rm = TRUE)

#scale the outcome
mctb_data$sooner_choice_scaled <- mctb_data$sooner_choice / sd_of_soonerchoice
mctb_data$endowment_later_scaled <- mctb_data$endowment_later / sd_of_soonerchoice

a_start <- 0.90
b_start <- 1
d_start <- 0.999

nls_ctb_model <- nls(sooner_choice ~ 
                       ( endowment_later * (b^t0 * d^k * pratio)^(1/(a-1)) ) / 
                       ( 1+pratio * (b^t0 * d^k * pratio)^(1/(a-1)) ),
          mctb_data, 
          start = list(a = a_start, 
                       b = b_start, 
                       d = d_start))

knitr::kable(broom::tidy(nls_ctb_model), digits = 4)
term estimate std.error statistic p.value
a 0.7451 0.0049 152.8705 0
b 0.9608 0.0045 213.8533 0
d 0.9986 0.0000 21015.3706 0
##STORE MODEL PARAMS##
alpha_nls <- summary(nls_ctb_model)$coefficients[1]
beta_nls<- summary(nls_ctb_model)$coefficients[2]
delta_nls <- summary(nls_ctb_model)$coefficients[3] 
mctb_data$alpha_nls <- alpha_nls
mctb_data$beta_nls <- beta_nls
mctb_data$delta_nls <- delta_nls

We can interpret these values as follows: \(a\) is considerably less than 1, indicating that generally people do not prefer the all-sooner or all-later options, but, rather, prefer to take some money sooner and some money later. The parameter \(b < 1\) indicates people discount even a bit more by a factor of \(b = 0.9608\) if they can get the sooner reward today, rather than waiting five weeks. In other words, for a constant delay of 5 weeks, they discount at a rate \(d=0.9986\) exponentiated by day then multiplied by \(0.9608\) if they get the sooner reward today (the sooner discount for a five week delay is \(0.9986^{5\times7} = 0.9523\), but if the sooner reward is today, it’s \(0.9608 \times 0.9986^{5\times7} = 0.915\)).

If we apply this to the full equation above, we can find the sooner payout that, combined with the portion of the payout received later, is equivalent to a $20 later payout for a participant with these parameter values (if the sooner payout is received today, for a price ratio of, say, 20/16 = 1.25):

round((sooner_choice_calc <- (20 * (beta_nls*delta_nls^(7*5)*1.25)^(1/(alpha_nls-1)) )/
  ( 1 + 1.25 * (beta_nls*delta_nls^(7*5)*1.25)^(1/(alpha_nls-1)) )),2)
round((amount_received_later <- 20-sooner_choice_calc*20/16), 2)
round(total_amount_received <- sooner_choice_calc + amount_received_later, 2)
## [1] 6.79
## [1] 11.51
## [1] 18.3

\[ 6.79 = \frac{20\cdot(0.9608^{1}\cdot 0.9986^{7\times5}\cdot1.25)^{\frac{1}{0.7451-1}}}{1+1.25\cdot(0.9608^{1}\cdot 0.9986^{7\times5}\cdot1.25)^{\frac{1}{0.7451-1}}} \]

So, these population parameters suggest that someone will be just as happy getting $6.79 today and $11.51 (for a total of $18.3) later as they would to receive the full $20 later.

Another way of expressing the discounting variable is in terms of a monthly discount rate:

(a_monthly_rate <- delta_nls^(-30) - 1)
## [1] 0.04275648

Individual-level parameters

Finding the cases that have no meaningful variability in their responses. We exclude trials where you could receive $20 sooner or $20 later.

odd_cases <- group_by(mctb_data, subject_id) %>%
  filter(soon_1 != late_6) %>%
  summarize(
    no_interior = case_when(
      any(c %in% 2:5) ~ FALSE,
      all(c %in% c(1,6)) ~ TRUE,
      TRUE ~ NA
    ),
    no_variance = case_when(
      all(c == 1) ~ TRUE,
      all(c == 6) ~ TRUE,
      TRUE ~ FALSE
    ))

print(paste0('Do all cases that have no variance also have no interior solutions? ', 
             all(odd_cases$no_interior[odd_cases$no_variance])))
## [1] "Do all cases that have no variance also have no interior solutions? TRUE"
#View(odd_cases)

table(odd_cases[,2:3])
##            no_variance
## no_interior FALSE TRUE
##       FALSE   533    0
##       TRUE    201  345

This indicates that 345 participants always answer either sooner or later exclusively (excluding the $20 now versus $20 later trials). 533 participants use the interior options, and 201 participants just switch between taking everything sooner, or everything later.

mctb_nlslist <- nlsList(
  sooner_choice_scaled ~ 
    ( endowment_later_scaled * (b^t0 * d^k_w * pratio)^(1/(a-1)) ) / 
    ( 1+pratio * (b^t0 * d^k_w * pratio)^(1/(a-1)) ) | subject_id, 
  data = select(mctb_data,
                t0, k_w, pratio, subject_id,
                sooner_choice_scaled, endowment_later_scaled), 
  start = c(a=alpha_nls,b=beta_nls, d=delta_nls),
  control = list(maxiter = 5000, minFactor = 1/2^30, tol = 1e-1))
## Warning: 492 errors caught in numericDeriv(form[[3L]], names(ind), env).  The error messages and their frequencies are
## 
## step factor 4.65661e-10 reduced below 'minFactor' of 9.31323e-10 
##                                                                9 
##                    number of iterations exceeded maximum of 5000 
##                                                               14 
##                                                singular gradient 
##                                                              100 
##  Missing value or an infinity produced when evaluating the model 
##                                                              369
converged <- lapply(mctb_nlslist,
                    function(amod) !is.null(amod))

converged_df <- data_frame(subject_id = as.numeric(names(converged)), 
                  nlslist_converged = unlist(converged),
                  nlslist = 1)

mctb_data_diagnostics <- left_join(
  left_join(mctb_data, odd_cases,
            by = 'subject_id'),
  converged_df[,-3],
  by = 'subject_id')

Take a look at responses from those with convergent versus non-convergent nls models. The first three panels are convergent, and the last three are non-convergnet.

ggplot2::theme_set(ggplot2::theme_minimal())
plot_mctb_responses(
  filter(mctb_data, subject_id %in% 
           c(converged_df$subject_id[converged_df$nlslist_converged][1:3],
             converged_df$subject_id[!converged_df$nlslist_converged][1:3])),
  adj.pratio = TRUE)

Use the nlslist object to start up the linear mixed effects model:

mctb_nlme.nlslist <- nlme.nlsList(model = mctb_nlslist,
                                  random = b + d ~ 1,
                                  groups = ~subject_id,
                                  method = 'ML')
summary(mctb_nlme.nlslist)
## Nonlinear mixed-effects model fit by maximum likelihood
##   Model: sooner_choice_scaled ~ (endowment_later_scaled * (b^t0 * d^k_w *      pratio)^(1/(a - 1)))/(1 + pratio * (b^t0 * d^k_w * pratio)^(1/(a -      1))) 
##  Data: select(mctb_data, t0, k_w, pratio, subject_id, sooner_choice_scaled,      endowment_later_scaled) 
##       AIC      BIC    logLik
##   51574.3 51631.44 -25780.15
## 
## Random effects:
##  Formula: list(b ~ 1, d ~ 1)
##  Level: subject_id
##  Structure: General positive-definite, Log-Cholesky parametrization
##          StdDev     Corr 
## b        0.13088065 b    
## d        0.02301261 0.704
## Residual 0.57988744      
## 
## Fixed effects: list(a ~ 1, b ~ 1, d ~ 1) 
##       Value   Std.Error    DF   t-value p-value
## a 0.9454517 0.000794510 24815 1189.9811       0
## b 0.9190270 0.004295194 24815  213.9663       0
## d 0.9829883 0.000715799 24815 1373.2737       0
##  Correlation: 
##   a      b     
## b -0.027       
## d -0.016  0.621
## 
## Standardized Within-Group Residuals:
##           Min            Q1           Med            Q3           Max 
## -4.578681e+00 -2.167673e-01 -1.581553e-05  2.787358e-01  4.000113e+00 
## 
## Number of Observations: 25896
## Number of Groups: 1079
qplot(sample = residuals(mctb_nlme.nlslist)/sd(residuals(mctb_nlme.nlslist)), 
      geom = c('qq', 'abline'), 
      main = "Q-Q plot for conditional residuals")

qplot(sample = ranef(mctb_nlme.nlslist)$b/sd(ranef(mctb_nlme.nlslist)$b), 
      geom = c('qq', 'abline'), 
      main = "Q-Q plot for  random b")

qplot(sample = ranef(mctb_nlme.nlslist)$d/sd(ranef(mctb_nlme.nlslist)$d), 
      geom = c('qq', 'abline'), 
      main = "Q-Q plot for  random d")

qplot(x = predict(mctb_nlme.nlslist), y = residuals(mctb_nlme.nlslist),
      geom = c('jitter', 'smooth'), width = .05, height = .05)
## Warning: Ignoring unknown parameters: height
## `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'

qplot(x = predict(mctb_nlme.nlslist), y = sqrt(abs(scale(residuals(mctb_nlme.nlslist)))),
      geom = c('jitter', 'smooth'), width = .05, height = .05)
## Warning: Ignoring unknown parameters: height
## `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'

plot(mctb_nlme.nlslist, subject_id ~ resid(.))

plot(mctb_nlme.nlslist, sooner_choice_scaled ~ fitted(.) | subject_id)
mctb_nlme.nlslist_coef_df <- as.data.frame(coef(mctb_nlme.nlslist))
mctb_nlme.nlslist_coef_df$subject_id <- as.numeric(rownames(mctb_nlme.nlslist_coef_df))

mctb_data_model_info <- left_join(mctb_data_diagnostics,
                                  mctb_nlme.nlslist_coef_df,
                                  by = 'subject_id')

mctb_model_results <- distinct(mctb_data_model_info,
                               subject_id, a, b, d, 
                               no_interior, no_variance, nlslist_converged)
write_csv(mctb_model_results, 
          path = file.path(data_dir_base, 
                           'JRP_CTB_model_paramsOSF.csv'))
mctb_pars_and_flags <- distinct(mctb_data_model_info, 
                  a, b, d, 
                  no_interior, no_variance, nlslist_converged,
                  subject_id)

flag_summaries <- group_by(mctb_pars_and_flags,
                           no_interior, no_variance) %>%
  summarize(a = mean(a, na.rm = T),
            b = mean(b, na.rm = T),
            d = mean(d, na.rm = T),
            n = n())

knitr::kable(bind_rows(flag_summaries,
                       as.data.frame(t(fixef(mctb_nlme.nlslist)))), digits = 4)
no_interior no_variance a b d n
FALSE FALSE 0.9455 0.9211 0.9838 533
TRUE FALSE 0.9455 0.9282 0.9851 201
TRUE TRUE 0.9455 0.9105 0.9805 345
NA NA 0.9455 0.9190 0.9830 NA
for(avar in list(list('b', .05), list('d', .005))) {
  print(ggplot(mctb_pars_and_flags,
         aes_string(x = avar[[1]], 
             group = 'nlslist_converged',
             fill = 'nlslist_converged')) + 
    geom_vline(xintercept = 1) +
    geom_histogram(position = position_dodge(),
                   aes(y = ..count..),
                   binwidth = avar[[2]]) + 
    facet_grid(no_interior ~ no_variance,
               labeller = labeller(
                 no_interior = as_labeller(c('TRUE' = 'No interior choices',
                                             'FALSE' = 'Has interior choices')),
                 no_variance = as_labeller(c('TRUE' = 'No choice variance',
                                             'FALSE' = 'Has choice variance'))),
               as.table = TRUE) + 
      labs(x = paste0('Parameter ', avar[[1]]),
           fill = 'Single-subject\nNLS model\nconverged?') + 
      theme_minimal())
}

Participants who have very low d values:

plot_small_d_sids <- unlist(distinct(filter(mctb_data_model_info, 
                                          d < .94, !no_variance),
                                   subject_id))
for(i in unique((1:23 -1) %/% 6)){
  last_index <- ifelse(i*6 < length(plot_small_d_sids),
                       (i+1)*6,
                       length(plot_small_d_sids))
  first_index <- i*6 + 1
  print(plot_mctb_responses(filter(
    mctb_data_model_info, 
    subject_id %in% plot_small_d_sids[first_index:last_index]))
    +theme_minimal())
}

Parameter-behavior correlation check

behavior_summary <- group_by(mctb_data, t0, subject_id) %>%
  summarize(sooner_choice = mean(sooner_choice)) %>%
  gather(key, value, -subject_id, -t0) %>%
  unite(key, key, t0) %>%
  spread(key, value) %>%
  mutate(mean_choice = (sooner_choice_0 + sooner_choice_1)/2,
         diff_choice = sooner_choice_1 - sooner_choice_0)

compare_df <- left_join(mctb_model_results, behavior_summary)
## Joining, by = "subject_id"
plot(compare_df$mean_choice, compare_df$d)

plot(compare_df$diff_choice, compare_df$b)

The above plots indicate a reasonable mapping of individual-level parameter values to the observed behavior.

Plotting all the choices

theme_set(theme_minimal())

mctb_sample_ids <- distinct(mctb_data, subject_id) %>%
  mutate(plot_group = (1:n()-1) %/% 5)

for(g in 0:max(mctb_sample_ids$plot_group)){
  mctb_data_sample <- left_join(
    filter(mctb_sample_ids, plot_group == g),
           mctb_data, by = 'subject_id')
  
  aplot <- plot_mctb_responses(mctb_data_sample)
  print(aplot)
}

Check sensitivity to letting \(\alpha\) vary over individuals

#this model shows poor convergence
mctb_nlme.nlslist_a <- nlme.nlsList(model = mctb_nlslist,
                                    random = a + b + d ~ 1,
                                    groups = ~subject_id,
                                    method = 'ML',
                                    control = nlmeControl(msVerbose = TRUE, 
                                                          niterEM = 100,
                                                          pnlsTol = 1.1),
                                    verbose = 2)
##   0:     119821.65:  1.72729  1.86088  3.52593 -1.24060 -12.2423 -25.7977
##   1:     119821.65:  1.72729  1.86088  3.52593 -1.24060 -12.2423 -25.7977
##   2:     119821.65:  1.72729  1.86088  3.52593 -1.24060 -12.2423 -25.7977
##   3:     119821.65:  1.72729  1.86088  3.52593 -1.24060 -12.2423 -25.7977
## 
## **Iteration 1
## LME step: Loglik: -24990.93, nlminb iterations: 3
## reStruct  parameters:
## subject_id1 subject_id2 subject_id3 subject_id4 subject_id5 subject_id6 
##    1.727291    1.860884    3.525932   -1.240604  -12.242275  -25.797727 
##  Beginning PNLS step: ..  completed fit_nlme() step.
## PNLS step: RSS =  5412.874 
##  fixed effects: 0.8425702  0.9113531  0.9825489  
##  iterations: 5 
## Convergence crit. (must all become <= tolerance = 1e-05):
##       fixed    reStruct 
##  0.03503821 12.64544046 
##   0:     115425.25: 0.380387  1.20965  2.81685 -0.914147 -0.897170 -12.8678
##   1:     115425.25: 0.380223  1.20944  2.81665 -0.914152 -0.897172 -12.8678
##   2:     115425.25: 0.380223  1.20944  2.81665 -0.914152 -0.897172 -12.8678
##   3:     115425.25: 0.380223  1.20944  2.81665 -0.914152 -0.897172 -12.8678
##   4:     115425.25: 0.380202  1.20948  2.81664 -0.914016 -0.897169 -12.8677
##   5:     115425.25: 0.380180  1.20942  2.81664 -0.913895 -0.897180 -12.8677
##   6:     115425.25: 0.380213  1.20944  2.81663 -0.913896 -0.897239 -12.8675
##   7:     115425.25: 0.380182  1.20946  2.81661 -0.913863 -0.897752 -12.8662
##   8:     115425.25: 0.380213  1.20936  2.81660 -0.913920 -0.898306 -12.8648
##   9:     115425.25: 0.380199  1.20934  2.81660 -0.913717 -0.898888 -12.8635
## 
## **Iteration 2
## LME step: Loglik: -20594.53, nlminb iterations: 9
## reStruct  parameters:
## subject_id1 subject_id2 subject_id3 subject_id4 subject_id5 subject_id6 
##   0.3801986   1.2093374   2.8165952  -0.9137169  -0.8988876 -12.8634672 
##  Beginning PNLS step: ..  completed fit_nlme() step.
## PNLS step: RSS =  5412.717 
##  fixed effects: 0.8425702  0.9113531  0.9825489  
##  iterations: 1 
## Convergence crit. (must all become <= tolerance = 1e-05):
##       fixed    reStruct 
## 0.000000000 0.002561383 
##   0:     115425.25: 0.380186  1.20932  2.81660 -0.913585 -0.899474 -12.8616
##   1:     115425.25: 0.380186  1.20932  2.81660 -0.913585 -0.899474 -12.8616
##   2:     115425.25: 0.380186  1.20932  2.81660 -0.913585 -0.899474 -12.8616
## 
## **Iteration 3
## LME step: Loglik: -20594.53, nlminb iterations: 2
## reStruct  parameters:
## subject_id1 subject_id2 subject_id3 subject_id4 subject_id5 subject_id6 
##   0.3801860   1.2093186   2.8165977  -0.9135847  -0.8994736 -12.8615728 
##  Beginning PNLS step: ..  completed fit_nlme() step.
## PNLS step: RSS =  5412.717 
##  fixed effects: 0.8425702  0.9113531  0.9825489  
##  iterations: 1 
## Convergence crit. (must all become <= tolerance = 1e-05):
##        fixed     reStruct 
## 0.000000e+00 1.495554e-07
summary(mctb_nlme.nlslist_a)
## Nonlinear mixed-effects model fit by maximum likelihood
##   Model: sooner_choice_scaled ~ (endowment_later_scaled * (b^t0 * d^k_w *      pratio)^(1/(a - 1)))/(1 + pratio * (b^t0 * d^k_w * pratio)^(1/(a -      1))) 
##  Data: select(mctb_data, t0, k_w, pratio, subject_id, sooner_choice_scaled,      endowment_later_scaled) 
##        AIC      BIC    logLik
##   41209.05 41290.67 -20594.53
## 
## Random effects:
##  Formula: list(a ~ 1, b ~ 1, d ~ 1)
##  Level: subject_id
##  Structure: General positive-definite, Log-Cholesky parametrization
##          StdDev     Corr       
## a        0.31186232 a     b    
## b        0.16056237 0.352      
## d        0.02550802 0.246 0.610
## Residual 0.42649076            
## 
## Fixed effects: list(a ~ 1, b ~ 1, d ~ 1) 
##       Value   Std.Error    DF   t-value p-value
## a 0.8425702 0.009861783 24815   85.4379       0
## b 0.9113531 0.005396787 24815  168.8696       0
## d 0.9825489 0.000824174 24815 1192.1624       0
##  Correlation: 
##   a     b    
## b 0.331      
## d 0.248 0.541
## 
## Standardized Within-Group Residuals:
##           Min            Q1           Med            Q3           Max 
## -5.644619e+00 -1.424034e-01 -7.529998e-20  4.850362e-01  5.602953e+00 
## 
## Number of Observations: 25896
## Number of Groups: 1079
summary(mctb_nlme.nlslist)
## Nonlinear mixed-effects model fit by maximum likelihood
##   Model: sooner_choice_scaled ~ (endowment_later_scaled * (b^t0 * d^k_w *      pratio)^(1/(a - 1)))/(1 + pratio * (b^t0 * d^k_w * pratio)^(1/(a -      1))) 
##  Data: select(mctb_data, t0, k_w, pratio, subject_id, sooner_choice_scaled,      endowment_later_scaled) 
##       AIC      BIC    logLik
##   51574.3 51631.44 -25780.15
## 
## Random effects:
##  Formula: list(b ~ 1, d ~ 1)
##  Level: subject_id
##  Structure: General positive-definite, Log-Cholesky parametrization
##          StdDev     Corr 
## b        0.13088065 b    
## d        0.02301261 0.704
## Residual 0.57988744      
## 
## Fixed effects: list(a ~ 1, b ~ 1, d ~ 1) 
##       Value   Std.Error    DF   t-value p-value
## a 0.9454517 0.000794510 24815 1189.9811       0
## b 0.9190270 0.004295194 24815  213.9663       0
## d 0.9829883 0.000715799 24815 1373.2737       0
##  Correlation: 
##   a      b     
## b -0.027       
## d -0.016  0.621
## 
## Standardized Within-Group Residuals:
##           Min            Q1           Med            Q3           Max 
## -4.578681e+00 -2.167673e-01 -1.581553e-05  2.787358e-01  4.000113e+00 
## 
## Number of Observations: 25896
## Number of Groups: 1079
plot(coef(mctb_nlme.nlslist_a)[,'b'], coef(mctb_nlme.nlslist)[,'b'])

plot(coef(mctb_nlme.nlslist_a)[,'d'], coef(mctb_nlme.nlslist)[,'d'])